/*******************************************************************************

This code file produces Figure 5, "Marginal Value of Public Funds of 421-a by Neighborhood Tabulation Area" and Figure A20, "Marginal Value of Public Funds of 421-a by Neighborhood Tabulation Area (High-WTP Scenario)."

Programs to run in HendrenSprungKeyser:
- metafile.do
- new_programs/run_program.ado
- new_programs/est_life_impact_nyc.ado
- new_programs/get_wtp_housing.ado
- ado/convert_rank_dollar.ado

For the high-WTP scenario, change line 162 in nyc421a.do.
For the close-moves scenario, change lines 17/18 in oi_tract_data.do

*******************************************************************************/

*** Manage settings

	run "$dir/code/modules/settings.do"
	
********************************************************************************
* Tempfile to store results
********************************************************************************	
	
	clear
	set obs 179
	gen nta = . 
	gen nta_code = ""
	gen NTAName = ""
	gen mvpf = .
	gen wtp = .
	gen mcost = .
	gen program_cost = .
	gen cost = .
	gen incidence = .
	
	tempfile mvpf_by_nta
	save `mvpf_by_nta', replace
	
********************************************************************************
* Obtain marginal cost per unit by NTA
********************************************************************************

	use "$data/clean/cleaned_data.dta", clear

	logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
	mat b = e(b)
	local coef_dtaxrate = b[1,1]
	
	predict xb, xb
	gen pr = exp(xb)/(1+exp(xb))
	gen pr_plus = exp(xb+`coef_dtaxrate'*0.001)/(1+exp(xb+`coef_dtaxrate'*0.001))
		
	gen assessprop = pr * (assesstot - assessland) * dtaxrate_onsite/sh_taxable
	gen assessprop_plus = pr_plus * (assesstot - assessland) * (dtaxrate_onsite/sh_taxable + 0.001)
	
	local sh_inclusionary = 0.20
	gen inclusionaryunits = `sh_inclusionary' * pr * unitsres 
	gen inclusionaryunits_plus = `sh_inclusionary' * pr_plus * unitsres
	
	collapse (sum) assessprop assessprop_plus inclusionaryunits inclusionaryunits_plus (mean) xcoord_latlong ycoord_latlong, by(nta_code nta)
		
	gen mcostperunit = (1-$beta) * (assessprop_plus-assessprop) / (inclusionaryunits_plus-inclusionaryunits)
	
	tempfile mcost
	save `mcost', replace
	
	use "$data/GIS_boundaries/neighborhood_boundaries/df.dta", clear
	merge 1:1 nta_code using `mcost'
	
	keep nta NTAName nta_code mcost
	
	save `mcost', replace
	
********************************************************************************
* Obtain 421-a incidence by NTA
********************************************************************************
	
	use "$data/clean/cleaned_data.dta", clear
	
	logit inclusionary_onsite dtaxrate_onsite $controls_lot $controls_block i.borough i.yearpermit, cl(nta)
		
	mat b = e(b)
	local coef_dtaxrate = b[1,1]
	predict xb, xb
	
	gen dlprofit = (1 / `coef_dtaxrate') * log(1 + exp(xb))
	gen mvalue = (assesstot-assessland)/(frac_assess*underassess)

	collapse (mean) dlprofit dtaxrate_onsite (rawsum) mvalue [aw=mval], by(borough nta nta_code)
	
	gen incidence = dlp/dtax
	
	replace incidence = 1 if incidence > 1 & !missing(incidence)
	
	bys borough: gegen inci_mean = mean(incidence) [aw=mvalue]
	replace incidence = inci_mean if missing(incidence)
	drop inci_mean borough
	
	keep nta_code incidence
	drop if missing(nta_code)
	
	tempfile incidence
	save `incidence', replace
	
********************************************************************************
* Merge together and run MVPF program
********************************************************************************
	
	preserve
	use "$data/raw/oi_pred_effect_by_nta_closemoves.dta", clear
	collapse (mean) pct_diff dist, by(nta_code)
	tempfile oi
	save `oi', replace
	restore
	
	merge 1:m nta_code using `oi', nogen
	merge m:1 nta_code using `mcost', nogen
	merge m:1 nta_code using `incidence', nogen
	
	keep if !missing(mcostperunit) & !missing(incidence) & !missing(pct_diff)
	
	local n = _N
	
	forvalues i = 1/`n' {
	
		di `i'
								
		global year_cost = mcost[`i']
		global pct_diff = pct_diff[`i']
		global incidence = incidence[`i']
		
		local ntaname = NTAName[`i']
		local nta = nta[`i']
		local nta_code = nta_code[`i']
		
		local mcost = mcost[`i']
		local incidence = incidence[`i']
	
		quietly run_program nyc421a
		local mvpf = $MVPF_nyc421a
		local wtp = $WTP_nyc421a
		local program_cost = $program_cost_nyc421a
		local cost = $cost_nyc421a
	
		preserve
		use `mvpf_by_nta', clear
	
		quietly replace mvpf = `mvpf' if _n == `i'
		quietly replace wtp = `wtp' if _n == `i'
		quietly replace program_cost = `program_cost' if _n == `i'
		quietly replace cost = `cost' if _n == `i'
		
		quietly replace mcost = `mcost' if _n == `i'
		quietly replace incidence = `incidence' if _n == `i'
		
		quietly replace NTAName = "`ntaname'" if _n == `i'
		quietly replace nta = `nta' if _n == `i'
		quietly replace nta_code = "`nta_code'" if _n == `i'
		
		quietly save `mvpf_by_nta', replace
		
		restore
	
	}
	
	use `mvpf_by_nta', clear
	drop if missing(nta_code)
		
	merge 1:1 nta_code using "$data/GIS_boundaries/neighborhood_boundaries/df.dta", nogen
	
	replace mvpf = 0 if wtp < 0
	replace mvpf = 9999 if mvpf < 0 & wtp > 0 & !missing(wtp)
	
	format mvpf %6.2f
			
	spmap mvpf using "$data/GIS_boundaries/neighborhood_boundaries/coords.dta", ///
		id(_ID) clm(custom) clb(0 0.5 1 5 10000) fcolor(Greys2) ndf(white) ///
		legend(order(5 "(5, {&infinity})" 4 "(1, 5]" 3 "(0.5, 1]" 2 "[0, 0.5]" 1 "No Data") position(4)) ///
		osize(vthin vthin vthin vthin vthin vthin) ndsize(vthin)
	
	graph export "$figs/mvpf_nta_lb_close.pdf", replace
	graph export "$figs_overleaf/mvpf_nta_lb_close.pdf", replace
	
